using Revise
using LatinHypercubeSampling
using DifferentialEquations
using Optim
using Statistics
using Distances
using Distributions
using Turing
using CSV
using DataFrames
using DataFramesMeta
using Chain: @chain
using Pipe: @pipe
using DataStructures
using Plots
using StatsPlots
using Interact
#pyplot()
Unable to load WebIO. Please make sure WebIO works for your Jupyter client. For troubleshooting, please see the WebIO/IJulia documentation.
includet("models/base.jl")
data = map(["wt", "ifnko"]) do geno
active = @linq CSV.read("data_old/data_$(geno)_active.csv", DataFrame) |>
transform(name="active") |>
select(genotype=:genotype, age=:age, name=:name, value=:active, weight=:weight)
counts = @linq CSV.read("data_old/data_$(geno)_counts.csv", DataFrame) |>
transform(name="total") |>
select(genotype=:genotype, age=:age, name=:name, value=:count, weight=:weight)
vcat([active, counts]...)
end
data = vcat(data...)
data = @linq data |> transform(link=[n .== "active" ? 2 : 1 for n in :name])
| genotype | age | name | value | weight | link | |
|---|---|---|---|---|---|---|
| String | Int64 | String | Float64 | Float64 | Int64 | |
| 1 | wt | 660 | active | 0.234043 | 676.37 | 2 |
| 2 | wt | 660 | active | 0.15 | 676.37 | 2 |
| 3 | wt | 660 | active | 0.2 | 676.37 | 2 |
| 4 | wt | 660 | active | 0.228571 | 676.37 | 2 |
| 5 | wt | 60 | active | 0.586207 | 42.1492 | 2 |
| 6 | wt | 60 | active | 0.467742 | 42.1492 | 2 |
| 7 | wt | 60 | active | 0.280702 | 42.1492 | 2 |
| 8 | wt | 60 | total | 764.0 | 1.97304e-5 | 1 |
| 9 | wt | 60 | total | 755.5 | 1.97304e-5 | 1 |
| 10 | wt | 60 | total | 1254.5 | 1.97304e-5 | 1 |
| 11 | wt | 60 | total | 845.5 | 1.97304e-5 | 1 |
| 12 | wt | 60 | total | 708.0 | 1.97304e-5 | 1 |
| 13 | wt | 60 | total | 1221.0 | 1.97304e-5 | 1 |
| 14 | wt | 60 | total | 804.0 | 1.97304e-5 | 1 |
| 15 | wt | 60 | total | 699.0 | 1.97304e-5 | 1 |
| 16 | wt | 210 | total | 232.0 | 0.000115616 | 1 |
| 17 | wt | 210 | total | 396.0 | 0.000115616 | 1 |
| 18 | wt | 210 | total | 238.0 | 0.000115616 | 1 |
| 19 | wt | 660 | total | 116.333 | 0.0001889 | 1 |
| 20 | wt | 660 | total | 88.3333 | 0.0001889 | 1 |
| 21 | wt | 660 | total | 226.0 | 0.0001889 | 1 |
| 22 | ifnko | 660 | active | 0.392857 | 646.154 | 2 |
| 23 | ifnko | 660 | active | 0.428571 | 646.154 | 2 |
| 24 | ifnko | 660 | active | 0.35 | 646.154 | 2 |
| 25 | ifnko | 60 | active | 0.333333 | 10227.9 | 2 |
| 26 | ifnko | 60 | active | 0.352941 | 10227.9 | 2 |
| 27 | ifnko | 60 | active | 0.340909 | 10227.9 | 2 |
| 28 | ifnko | 60 | total | 584.0 | 2.12442e-5 | 1 |
| 29 | ifnko | 60 | total | 691.0 | 2.12442e-5 | 1 |
| 30 | ifnko | 60 | total | 1022.5 | 2.12442e-5 | 1 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
sim_ifnko_both = CSV.read("data_old/solution_ifnko_self_renewal_and_quiescence.csv", DataFrame)
sim_wt_both = CSV.read("data_old/solution_wt_self_renewal_and_quiescence.csv", DataFrame)
sims = vcat(@transform(sim_wt_both, genotype="wt"), @transform(sim_ifnko_both, genotype="ko"))
| t | active | counts | b | r | genotype | |
|---|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Float64 | String | |
| 1 | 0.0 | 0.474127 | 1986.31 | 0.482161 | 0.842565 | wt |
| 2 | 0.157915 | 0.474119 | 1981.27 | 0.482174 | 0.842313 | wt |
| 3 | 0.271147 | 0.474101 | 1977.67 | 0.482184 | 0.842133 | wt |
| 4 | 0.423967 | 0.474068 | 1972.82 | 0.482197 | 0.84189 | wt |
| 5 | 0.599921 | 0.474017 | 1967.26 | 0.482212 | 0.84161 | wt |
| 6 | 0.804687 | 0.473946 | 1960.82 | 0.48223 | 0.841284 | wt |
| 7 | 1.0378 | 0.473854 | 1953.52 | 0.48225 | 0.840913 | wt |
| 8 | 1.29979 | 0.473741 | 1945.36 | 0.482273 | 0.840497 | wt |
| 9 | 1.59052 | 0.473609 | 1936.35 | 0.482298 | 0.840035 | wt |
| 10 | 1.90964 | 0.473459 | 1926.54 | 0.482325 | 0.839528 | wt |
| 11 | 2.25654 | 0.473291 | 1915.94 | 0.482355 | 0.838978 | wt |
| 12 | 2.63039 | 0.473109 | 1904.62 | 0.482387 | 0.838385 | wt |
| 13 | 3.03023 | 0.472912 | 1892.6 | 0.482421 | 0.837752 | wt |
| 14 | 3.45497 | 0.472702 | 1879.95 | 0.482458 | 0.837079 | wt |
| 15 | 3.90343 | 0.47248 | 1866.72 | 0.482496 | 0.83637 | wt |
| 16 | 4.37437 | 0.472247 | 1852.97 | 0.482536 | 0.835626 | wt |
| 17 | 4.86656 | 0.472002 | 1838.74 | 0.482578 | 0.834848 | wt |
| 18 | 5.37872 | 0.471748 | 1824.09 | 0.482621 | 0.83404 | wt |
| 19 | 5.9096 | 0.471485 | 1809.08 | 0.482666 | 0.833204 | wt |
| 20 | 6.458 | 0.471213 | 1793.75 | 0.482712 | 0.83234 | wt |
| 21 | 7.02271 | 0.470933 | 1778.15 | 0.482759 | 0.831452 | wt |
| 22 | 7.6026 | 0.470646 | 1762.33 | 0.482808 | 0.830541 | wt |
| 23 | 8.19658 | 0.470351 | 1746.32 | 0.482857 | 0.829609 | wt |
| 24 | 8.80361 | 0.470051 | 1730.17 | 0.482908 | 0.828658 | wt |
| 25 | 9.42271 | 0.469744 | 1713.91 | 0.482959 | 0.827688 | wt |
| 26 | 10.053 | 0.469432 | 1697.58 | 0.483011 | 0.826703 | wt |
| 27 | 10.6935 | 0.469115 | 1681.19 | 0.483064 | 0.825702 | wt |
| 28 | 11.3435 | 0.468793 | 1664.8 | 0.483118 | 0.824688 | wt |
| 29 | 12.0023 | 0.468467 | 1648.4 | 0.483172 | 0.823662 | wt |
| 30 | 12.669 | 0.468138 | 1632.04 | 0.483226 | 0.822624 | wt |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
I want a model that shares one $b(\text{NSC})$ for both populations, and I want to see if that fits.
For workability I will use a $b(\text{NSC})$ instead of the other canditates. This will make it much easier to estimate initial values (as $b$ is dependent on $\text{NSC}_0$ and not the ratio which we need to estimate)
p = plot(xlab="NSC", ylab="b", xflip=false)
for group in groupby(sims, :genotype)
@df group plot!(p, :counts, :b, label=first(group.genotype))
end
p
What to choose as a shape for $b$?
Negative log10 has some good support, as this seems to be the relationship in these two fits:
p = plot(xlab="NSC", ylab="b", xflip=false, xscale=:log10)
for group in groupby(sims, :genotype)
@df group plot!(p, :counts, :b, label=first(group.genotype))
end
p
This does however have the problem of being technically able to go $ b>0.5 $, which would. In reality I do think that it will however stop before then, as the stem cell population has sufficiently decreased.
counts = @where(sims, :genotype .== "ko").counts
b = @where(sims, :genotype .== "ko").b
opt = optimize(x -> euclidean(x[1].+x[2]*log10.(counts), b), [1.0, 1.0])
* Status: success
* Candidate solution
Final objective value: 7.017176e-04
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 53
f(x) calls: 104
opt.minimizer
2-element Vector{Float64}:
0.5985516865001772
-0.045088147649260335
p = plot(xlab="NSC", ylab="b", xflip=false, xscale=:log10)
for group in groupby(sims, :genotype)
@df group plot!(p, :counts, :b, label=first(group.genotype), lw=2)
end
plot!(x -> opt.minimizer[1]+(opt.minimizer[2]*log10(x)), lab="ex", ls=:dash, lw=2)
p
Unfortunately this log-linear relationship can produce negative values of $b$ given sufficient $\text{NSC}$. As $b$ is a probability, this is not really meaningful and thus this log-linear relationship is only really viable as an approximation with sufficiently small $\text{NSC}$ or a suficient set of parameters.
A sigmoid is another viable alternative for this relationship. In the simulations we have thus far we have only ever explored one half of this sigmoid and thus there isn't really too much evidence for this in fact being sigmoidally shaped. A sigmoid is however supported by the argument that $b<0$ is meaningless and thus we need a function that (similar to the previous function) $$ [0,\infty) \rightarrow [0, 0.5] $$ or something close $$ [0,\infty) \rightarrow [0, 0.5) $$
The most convenient example of this is probably the Hill function (as this has use in biochemical modelling and is frequently used with some concentration as a support and a fraction as an output, in essence the same type of function we are looking for here).
It appears to be expressed as: $$ y = \frac{x^n}{K_d+x^n} $$ With ligand concentrations $x$, fraction bound $y$, dissociation constant $K_d$ and Hill coefficient $n$. This gives us two parameters to determine the shape: $K_d$ determines the general location of the function and $n$ determines the 'steepness' of it.
hill(kd, n, x) = (x^n)/(kd+x^n)
hill(kd, n) = x -> hill(kd, n, x)
hill (generic function with 2 methods)
This second parameterisation is nicer as $K_a$ represents the half-occupancy concentration, or for us the number of stem cells needed for $b=0.5$, i.e. the steady state or the number of stem cells we will plateau at.
hill(ka, n, x) = 1/(1+(ka/x)^n)
hill(ka, n) = x -> hill(ka, n, x)
hill (generic function with 2 methods)
points = 10.0 .^ (-3:0.1:5)
plot(xscale=:log10, xlim=(0.001, 1000))
plot!(points, hill(.001, 1), lab="")
plot!(points, hill(.01, 1), lab="")
plot!(points, hill(.1, 1), lab="")
plot!(points, hill(1, 1), lab="")
plot!(points, hill(10, 1), lab="")
plot!(points, hill(100, 1), lab="")
plot!(points, hill(1000, 1), lab="")
points = 10.0 .^ (-3:0.1:5)
plot(xscale=:log10, xlim=(0.001, 1000))
plot!(points, hill(1, 1/4), lab="")
plot!(points, hill(1, 1/2), lab="")
plot!(points, hill(1, 1), lab="")
plot!(points, hill(1, 2), lab="")
plot!(points, hill(1, 4), lab="")
There is two possible problems with this:
p = plot(xlab="NSC", ylab="b", xflip=false)
for group in groupby(sims, :genotype)
@df group plot!(p, :counts, :b, label=first(group.genotype))
end
p
We can of course easily fix the second issue:
plot(points, x -> 1-hill(1, 1, x), lab="", xscale=:log10)
And eye-balling a fit, it seems to be plausible that we can fit this function:
p = plot(points, x -> 1-hill(150, 0.05, x), lab="", xscale=:log10)
for group in groupby(sims, :genotype)
@df group plot!(p, :counts, :b, label=first(group.genotype))
end
p
Although I would really like this to be parameterised with a slope instead, this will do for now.
counts = @where(sims, :genotype .== "ko").counts
b = @where(sims, :genotype .== "ko").b
koopt = optimize(x -> euclidean(1 .- hill.(x[1], x[2], counts), b), [1.0, 1.0])
* Status: success
* Candidate solution
Final objective value: 3.478333e-04
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 61
f(x) calls: 120
counts = @where(sims, :genotype .== "wt").counts
b = @where(sims, :genotype .== "wt").b
wtopt = optimize(x -> euclidean(1 .- hill.(x[1], x[2], counts), b), [1.0, 1.0])
* Status: success
* Candidate solution
Final objective value: 1.184046e-02
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 57
f(x) calls: 109
hcat(koopt.minimizer, wtopt.minimizer)
2×2 Matrix{Float64}:
153.555 119.102
0.0785792 0.0267004
points = 10.0 .^ (-1:0.1:4)
p = plot(xscale=:log10,)
plot!(p, points, x -> 1-hill(koopt.minimizer[1], koopt.minimizer[2], x), lab="", lc=:gray, ls=:dash)
plot!(p, points, x -> 1-hill(wtopt.minimizer[1], wtopt.minimizer[2], x), lab="", lc=:gray, ls=:dash)
for (i, group) in enumerate(groupby(sims, :genotype))
@df group plot!(p, :counts, :b, label=first(group.genotype), lw=2, lc=i)
end
p
points = 10.0 .^ (2:0.1:3.5)
p = plot(xscale=:log10,)
plot!(p, points, x -> 1-hill(koopt.minimizer[1], koopt.minimizer[2], x), lab="", lc=:gray, ls=:dash)
plot!(p, points, x -> 1-hill(wtopt.minimizer[1], wtopt.minimizer[2], x), lab="", lc=:gray, ls=:dash)
for (i, group) in enumerate(groupby(sims, :genotype))
@df group plot!(p, :counts, :b, label=first(group.genotype), lw=2, lc=i)
end
p
struct PopSelf <: ODEModel
tspan::Tuple{Float64, Float64}
times::Vector{Float64}
values::Vector{Int64}
fixed::Vector{Pair{Int64, Float64}} # pointer to thing to fix and value
PopSelf(tspan, times, values,
fixed::Vector{Pair{I, N}}) where {I <: Integer, N <: Number} = new(tspan, times, values, fixed)
end
function ratefun(model::PopSelf)
function(du, u, p, t)
_, r₀, βr, ba, bb, pₛ = p
Q = u[1]; A = u[2]
b = 1-hill(ba, bb, Q+A)
du[1] = dQ = -exponential_decay(r₀, βr, t) * Q + 2 * b * pₛ * A
du[2] = dA = exponential_decay(r₀, βr, t) * Q - pₛ * A
end
end
function initial(t::PopSelf, x::AbstractVector)
nsc₀, r₀, _, ba, bb, pₛ = x
b₀ = 1-hill(ba, bb, nsc₀)
ratio = sqrt(((pₛ - r₀)/(2*r₀))^2 + (2*b₀*pₛ) / r₀)
nsc₀ .* [1-1/(ratio+1), 1/(ratio+1)]
end
link(t::PopSelf, x::AbstractArray) = hcat(x[:,1] .+ x[:,2], x[:,2] ./ (x[:,1] .+ x[:,2]))
parameter_names(t::Type{PopSelf}) = [:nsc₀, :r₀, :βr, :ba, :bb, :pₛ]
bounds(t::PopSelf) = [(100.0, 3000.0), (0.0, 1.0), (0.0, 0.1), (10.0, 1000.0), (0.005, 2.0), (0.0, 1.0)]
output_names(t::PopSelf) = ["qNSC", "aNSC"]
link_names(t::PopSelf) = ["Total NSC", "Fraction active NSC"]
link_names (generic function with 5 methods)
data_wt = @where(data, :genotype .== "wt")
model_wt = PopSelf((0.0, 700.0), Float64.(data_wt.age), data_wt.link, Dict(:pₛ => 0.95))
PopSelf((0.0, 700.0), [660.0, 660.0, 660.0, 660.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0 … 60.0, 60.0, 60.0, 60.0, 210.0, 210.0, 210.0, 660.0, 660.0, 660.0], [2, 2, 2, 2, 2, 2, 2, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [6 => 0.95])
model_starts = starts(model_wt, n=28*5)
140×5 Matrix{Float64}:
558.993 0.057554 0.0230216 736.475 1.44025
1205.76 0.0215827 0.0107914 287.77 0.808741
2582.73 0.294964 0.0489209 152.446 1.68424
1184.89 0.0863309 0.0244604 750.719 1.46896
204.317 0.280576 0.0604317 323.381 0.708273
1518.71 0.820144 0.0532374 572.662 0.320755
2207.19 0.129496 0.0410072 66.9784 0.665216
1873.38 0.625899 0.00215827 757.842 1.05273
1268.35 0.76259 0.0215827 636.763 1.02403
2979.14 0.640288 0.0345324 601.151 0.679568
1247.48 0.366906 0.0258993 394.604 0.0767626
934.532 0.546763 0.094964 729.353 1.41155
1977.7 0.215827 0.018705 95.4676 1.36849
⋮
1059.71 0.877698 0.0294964 123.957 1.66989
2812.23 0.611511 0.00503597 501.439 0.134173
997.122 0.985612 0.0582734 408.849 0.894856
2061.15 0.992806 0.0669065 857.554 0.363813
2290.65 0.338129 0.0561151 800.576 1.08144
1476.98 0.733813 0.00719424 230.791 0.34946
266.906 0.374101 0.0920863 957.266 0.995324
2311.51 0.661871 0.0884892 444.46 1.98565
2561.87 0.748201 0.0625899 864.676 1.51201
684.173 0.316547 0.0482014 764.964 1.29673
2603.6 0.532374 0.0330935 793.453 1.62683
1581.29 0.352518 0.0309353 337.626 0.794388
model_wt_opt = optimise(model_wt, (x, y) -> weuclidean(x, y, data_wt.weight), data_wt.value, model_starts);
minimae = minimum.(model_wt_opt)
min_index = argmin(minimae)
#scatter(sort(minimae), lab="", mc=:black)
scatter(minimae, lab="", mc=:black)
scatter!([min_index], [minimum(model_wt_opt[min_index])], lab="", mc=:red)
best_wt_opt = model_wt_opt[argmin(minimum.(model_wt_opt))]
* Status: success
* Candidate solution
Final objective value: 4.001431e+00
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 2 (vs limit Inf)
Iterations: 793
f(x) calls: 1361
minimum(best_wt_opt)
4.001430878361646
model_wt_plot = dense(model_wt)
param_wt = parameter_dict(model_wt,transform(model_wt, best_wt_opt.minimizer))
Dict{Symbol, Float64} with 6 entries:
:βr => 0.00182219
:r₀ => 0.805172
:nsc₀ => 2169.25
:bb => 0.0295317
:ba => 124.3
:pₛ => 0.95
plotlyjs()
p1 = plot(yscale=:log10, ylab="Total NSC", xlab="Time (days)")
p2 = plot(ylab="Fraction active", xlab="Time (days)")
for opt in filter(x -> minimum(x) < 7.5, model_wt_opt)
sim = parameter_array(model_wt_plot, transform(model_wt_plot, opt.minimizer)) |> simulate(model_wt_plot) |> link(model_wt_plot)
plot!(p1, model_wt_plot.times, sim[:,1], lab="", lc=:gray, lw=1, alpha=0.4)
plot!(p2, model_wt_plot.times, sim[:,2], lab="", lc=:gray, lw=1, alpha=0.4)
end
#best_i = 2
#sim = parameter_array(model_wt_plot, transform(model_wt_plot, model_wt_opt[best_i].minimizer)) |> simulate(model_wt_plot) |> link(model_wt_plot)
#plot!(p1, model_wt_plot.times, sim[:,1], lab="", lc=:green, lw=2)
#plot!(p2, model_wt_plot.times, sim[:,2], lab="", lc=:green, lw=2)
sim = param_wt |> simulate(model_wt_plot) |> link(model_wt_plot)
plot!(p1, model_wt_plot.times, sim[:,1], lab="", lc=:red, lw=2)
plot!(p2, model_wt_plot.times, sim[:,2], lab="", lc=:red, lw=2)
@df @where(data_wt, :link .== 1) scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data_wt, :link .== 2) scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2)
ba = param_wt[:ba]
bb = param_wt[:bb]
##nsc₀, r₀, _, ba, bb, pₛ = x
b = 1 .- hill.(ba, bb, sim[:,1])
p = plot(xflip=false, legend=:topright)#, xscale=:log10)
plot!(p, sim[:,1], b, lab="new", lw=2)
for (i, group) in enumerate(groupby(sims, :genotype))
@df group plot!(p, :counts, :b, label=first(group.genotype), lw=2)
end
p
data_ko = @where(data, :genotype .== "ifnko")
model_ko = PopSelf((0.0, 700.0), Float64.(data_ko.age), data_ko.link, Dict(:pₛ => 0.95))
PopSelf((0.0, 700.0), [660.0, 660.0, 660.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 180.0, 180.0, 180.0, 660.0, 660.0, 660.0], [2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [6 => 0.95])
model_starts = starts(model_ko, n=28*5)
140×5 Matrix{Float64}:
2394.96 0.841727 0.0460432 437.338 1.55507
725.899 0.0863309 0.042446 508.561 0.29205
1143.17 0.928058 0.0956835 202.302 1.48331
1727.34 0.956835 0.0431655 230.791 1.81342
2186.33 0.611511 0.0158273 88.3453 0.320755
2687.05 0.00719424 0.0920863 494.317 1.03838
2207.19 0.136691 0.0669065 829.065 1.64119
1164.03 0.892086 0.0834532 636.763 0.0480576
183.453 0.848921 0.0510791 180.935 1.39719
2728.78 0.647482 0.0928058 401.727 0.966619
538.129 0.676259 0.0791367 964.388 1.31108
350.36 0.165468 0.0863309 423.094 0.937914
1080.58 0.913669 0.0848921 764.964 1.7273
⋮
2228.06 0.0215827 0.0654676 572.662 0.507338
1351.8 0.230216 0.0561151 245.036 1.13885
1873.38 0.985612 0.0942446 836.187 0.363813
2603.6 0.345324 0.0345324 921.655 0.492986
2582.73 0.42446 0.0165468 24.2446 0.952266
788.489 0.877698 0.0992806 273.525 0.34946
2770.5 0.201439 0.094964 707.986 0.808741
517.266 0.158273 0.00863309 743.597 0.607806
412.95 0.381295 0.0395683 992.878 1.4546
663.309 0.76259 0.0827338 66.9784 0.923561
2290.65 0.215827 0.0705036 302.014 1.94259
2520.14 0.935252 0.0884892 472.95 0.880504
model_ko_opt = optimise(model_ko, (x, y) -> weuclidean(x, y, data_ko.weight), data_ko.value, model_starts);
minimae = minimum.(model_ko_opt)
min_index = argmin(minimae)
#scatter(sort(minimae), lab="", mc=:black)
scatter(minimae, lab="", mc=:black)
scatter!([min_index], [minimum(model_ko_opt[min_index])], lab="", mc=:red)
best_ko_opt = model_ko_opt[argmin(minimum.(model_ko_opt))]
* Status: success
* Candidate solution
Final objective value: 4.494265e+00
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 385
f(x) calls: 660
minimum(best_ko_opt)
4.494265394470573
model_ko_plot = dense(model_ko)
param_ko = parameter_dict(model_ko, transform(model_ko, best_ko_opt.minimizer))
Dict{Symbol, Float64} with 6 entries:
:βr => 1.51822e-21
:r₀ => 0.492025
:nsc₀ => 3000.0
:bb => 0.0847458
:ba => 155.637
:pₛ => 0.95
plotlyjs()
p1 = plot(yscale=:log10, ylab="Total NSC", xlab="Time (days)")
p2 = plot(ylab="Fraction active", xlab="Time (days)")
for opt in filter(x -> minimum(x) < 7.5, model_ko_opt)
sim = parameter_array(model_ko_plot, transform(model_ko_plot, opt.minimizer)) |> simulate(model_ko_plot) |> link(model_ko_plot)
plot!(p1, model_ko_plot.times, sim[:,1], lab="", lc=:gray, lw=1, alpha=0.4)
plot!(p2, model_ko_plot.times, sim[:,2], lab="", lc=:gray, lw=1, alpha=0.4)
end
#best_i = 2
#sim = parameter_array(model_ko_plot, transform(model_ko_plot, model_ko_opt[best_i].minimizer)) |> simulate(model_ko_plot) |> link(model_ko_plot)
#plot!(p1, model_ko_plot.times, sim[:,1], lab="", lc=:green, lw=2)
#plot!(p2, model_ko_plot.times, sim[:,2], lab="", lc=:green, lw=2)
sim = param_ko |> simulate(model_ko_plot) |> link(model_ko_plot)
plot!(p1, model_ko_plot.times, sim[:,1], lab="", lc=:red, lw=2)
plot!(p2, model_ko_plot.times, sim[:,2], lab="", lc=:red, lw=2)
@df @where(data_ko, :link .== 1) scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data_ko, :link .== 2) scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2)
ba = param_ko[:ba]
bb = param_ko[:bb]
##nsc₀, r₀, _, ba, bb, pₛ = x
b = 1 .- hill.(ba, bb, sim[:,1])
p = plot(xflip=false, legend=:topright)#, xscale=:log10)
plot!(p, sim[:,1], b, lab="new", lw=2)
for (i, group) in enumerate(groupby(sims, :genotype))
@df group plot!(p, :counts, :b, label=first(group.genotype), lw=2)
end
p
ba = param_ko[:ba]
bb = param_ko[:bb]
b_ko = 1 .- hill.(ba, bb, sim[:,1])
ba = param_wt[:ba]
bb = param_wt[:bb]
b_wt = 1 .- hill.(ba, bb, sim[:,1])
p = plot(xflip=false, legend=:topright)#, xscale=:log10)
plot!(p, sim[:,1], b_ko, lab="ko new", lw=2)
plot!(p, sim[:,1], b_wt, lab="wt new", lw=2)
for (i, group) in enumerate(groupby(sims, :genotype))
@df group plot!(p, :counts, :b, label=first(group.genotype), lw=2)
end
p
ba = param_ko[:ba]
bb = param_ko[:bb]
b_ko = 1 .- hill.(ba, bb, sim[:,1])
ba = param_wt[:ba]
bb = param_wt[:bb]
b_wt = 1 .- hill.(ba, bb, sim[:,1])
p = plot(xflip=false, legend=:topright)#, xscale=:log10)
plot!(p, model_ko_plot.times, b_ko, lab="ko new", lw=2)
plot!(p, model_wt_plot.times, b_wt, lab="wt new", lw=2)
for (i, group) in enumerate(groupby(sims, :genotype))
@df group plot!(p, :t, :b, label=first(group.genotype), lw=2)
end
p
param_ko
Dict{Symbol, Float64} with 6 entries:
:βr => 1.27258e-13
:r₀ => 0.492023
:nsc₀ => 3000.0
:bb => 0.0847444
:ba => 155.637
:pₛ => 0.95
param_fix = Dict(
:pₛ => 0.95,
:ba => param_ko[:ba],
:bb => param_ko[:bb],
)
model_wt_fixb = PopSelf((0.0, 700.0), Float64.(data_wt.age), data_wt.link, param_fix)
PopSelf((0.0, 700.0), [660.0, 660.0, 660.0, 660.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0 … 60.0, 60.0, 60.0, 60.0, 210.0, 210.0, 210.0, 660.0, 660.0, 660.0], [2, 2, 2, 2, 2, 2, 2, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [5 => 0.08474442921555739, 4 => 155.6366817708986, 6 => 0.95])
model_starts = starts(model_wt_fixb, n=28*5)
140×3 Matrix{Float64}:
1873.38 0.920863 0.0669065
1351.8 0.316547 0.0395683
308.633 0.52518 0.0539568
1310.07 0.705036 0.0697842
2228.06 1.0e-10 0.0316547
2311.51 0.755396 0.0769784
1393.53 0.94964 0.0136691
1935.97 0.309353 0.0158273
1789.93 0.0935252 0.0230216
3000.0 0.899281 0.0446043
830.216 0.280576 0.0755396
1810.79 0.618705 0.023741
1998.56 0.503597 0.0884892
⋮
475.54 0.496403 0.0251799
454.676 0.654676 0.00143885
162.59 0.834532 0.0280576
809.353 0.935252 0.028777
976.259 0.719424 0.0913669
2937.41 0.230216 0.057554
1664.75 0.208633 0.0942446
1247.48 0.345324 0.0848921
997.122 0.985612 0.0165468
2165.47 0.848921 0.0366906
851.079 0.172662 0.00791367
1205.76 0.669065 0.0201439
model_wt_fixb_opt = optimise(model_wt_fixb, (x, y) -> weuclidean(x, y, data_wt.weight), data_wt.value, model_starts);
best_wt_fixb_opt = model_wt_fixb_opt[argmin(minimum.(model_wt_fixb_opt))]
* Status: success
* Candidate solution
Final objective value: 4.795049e+00
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 73
f(x) calls: 144
minimum(best_wt_fixb_opt)
4.7950493837089585
model_wt_fixb_plot = dense(model_wt_fixb)
param_fixb_wt = parameter_dict(model_wt_fixb, transform(model_wt_fixb, best_wt_fixb_opt.minimizer))
Dict{Symbol, Float64} with 6 entries:
:βr => 0.000195959
:r₀ => 0.275258
:nsc₀ => 3000.0
:bb => 0.0847444
:ba => 155.637
:pₛ => 0.95
plotlyjs()
p1 = plot(yscale=:log10, ylab="Total NSC", xlab="Time (days)")
p2 = plot(ylab="Fraction active", xlab="Time (days)")
for opt in model_wt_fixb_opt #filter(x -> minimum(x) < 7.5, model_wt_fixb_opt)
sim = parameter_array(model_wt_fixb_plot, transform(model_wt_fixb_plot, opt.minimizer)) |> simulate(model_wt_fixb_plot) |> link(model_wt_fixb_plot)
plot!(p1, model_wt_fixb_plot.times, sim[:,1], lab="", lc=:gray, lw=1, alpha=0.4)
plot!(p2, model_wt_fixb_plot.times, sim[:,2], lab="", lc=:gray, lw=1, alpha=0.4)
end
best_i = 2
sim = parameter_array(model_wt_fixb_plot, transform(model_wt_fixb_plot, model_wt_opt[best_i].minimizer)) |> simulate(model_wt_fixb_plot) |> link(model_wt_fixb_plot)
plot!(p1, model_wt_fixb_plot.times, sim[:,1], lab="", lc=:green, lw=2)
plot!(p2, model_wt_fixb_plot.times, sim[:,2], lab="", lc=:green, lw=2)
sim = param_fixb_wt |> simulate(model_wt_fixb_plot) |> link(model_wt_fixb_plot)
plot!(p1, model_wt_fixb_plot.times, sim[:,1], lab="", lc=:red, lw=2)
plot!(p2, model_wt_fixb_plot.times, sim[:,2], lab="", lc=:red, lw=2)
@df @where(data_wt, :link .== 1) scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data_wt, :link .== 2) scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2)